!     written by Paul Colucci
!     CFD Laboratory
!     University at Buffalo
module params
  implicit none
  save
  
  integer, parameter :: nx = 256
  integer, parameter :: ny = 128
  integer, parameter :: ns = 3
  integer, parameter :: neq = 4 + (ns-1)
  integer, parameter :: nse = ns-1

end module params

module var
  use params
  implicit none
  save
  
  real, dimension(nx,ny) :: rho, u, v, p, T, E

end module var

module strs
  use params
  implicit none
  save
  
  real, dimension(nx,ny) :: Txx, Txy, Tyy, qx, qy

end module strs

module species
  use params
  implicit none
  save
  
  real, dimension(nx,ny,ns)  :: fs
  real, dimension(nx,ny,nse) :: mfluxx, mfluxy

end module species

module vec
  use params
  implicit none
  save
  
  real, dimension(nx,ny,neq) :: A, As, Aold
  real, dimension(nx,ny)     :: F, G, H

end module vec

module coord
  use params
  implicit none
  save
  
  real, dimension(nx,ny) :: x, y

end module coord

module trnsfrm
  use params
  implicit none
  save
  
  real                     :: dxi, deta
  real, dimension(nx,ny)   :: rjaci
  real, dimension(nx,ny,4) :: rmet, rmeti

end module trnsfrm

module steps
  implicit none
  save
  
  real    :: time, dt, stable
  integer :: iter

end module steps

module props
  implicit none
  save
  
  real :: mu, gamma, Pr, Sc, Le

end module props

module react
  use params
  implicit none
  save
  
  real                   :: Da1, Ce, Ze
  real, dimension(nx,ny) :: omega

end module react

module same
  implicit none
  save
  
  real :: Re1, M1

end module same

module delta
  use params
  implicit none
  save
  
  real, dimension(nx,ny) :: dx, dy

end module delta

module veccap
  use params
  implicit none
  save

  real, dimension(nx,ny) :: Fc, Gc

end module veccap
!=====================

program colucci
  use params
  use var
  use strs
  use species
  use vec
  use coord
  use trnsfrm
  use steps
  use props
  use react
  use same
  implicit none

  integer :: itmax, incrmnt, i, j, k, nframes 

  open(20,file='grid.dat') 
  open(21,file='soln.dat') 
  open(18,file='movie.dat') 
  open(19,file='movie.xyz') 

  itmax=4000    
  gamma=1.4  
  M1=2.4  
  Re1=400.  
  mu=1.    
  Pr=1. 
  Sc=1. 
  Le=1./Pr/Sc 
  Da1=5.    
  Ce=0.0  
  Ze=0.0 
  stable=0.5    
  !     stable is max CFL and is scheme dependent
  !     MacCormack           max cfl=1.0           typical=0.9
  !     GT                   max cfl=2./3          typical=0.6
  !     DCPS (compact)       max cfl=1./sqrt(3.)   typical=0.5

  !     Read external grid and compute metrics numerically 
  call rgrid 
  call metric 

  !     Initial conditions  
  call init 

  time=0.0 
  iter=0
  !     temporal integration of governing equations 
  nframes=10 
  incrmnt=itmax/nframes    
  write(18,*) nx/4,ny/4,nframes,1
1000 iter=iter+1 
  call tstep 
  time=time+dt     

  if (mod(iter,100).eq.0) write(*,*) iter,time,dt  

  call advance 

  !-----write product to file at desired timesteps for movie

  if (mod(iter,incrmnt).eq.0) then
     write(18,*)((fs(i,j,3),i=1,nx,4),j=1,ny,4)
  endif
  !     end of time loop      
  if (iter.lt.itmax) goto 1000 

  !     write grid to output file

  write(20,*) nx/2,ny/2,1  
  write(20,*) ((x(i,j),i=1,nx,2),j=1,ny,2), &
       ((y(i,j),i=1,nx,2),j=1,ny,2), &
       (( 0. ,i=1,nx,2),j=1,ny,2) 

  write(19,*) nx/4,ny/4,nframes 
  write(19,*) (((x(i,j),i=1,nx,4),j=1,ny,4),k=1,nframes), &
       (((y(i,j),i=1,nx,4),j=1,ny,4),k=1,nframes), &
       ((( 0. ,i=1,nx,4),j=1,ny,4),k=1,nframes)      
  !     write data to output file 

  write(21,*) nx/2,ny/2,1 
  write(21,*) 0,0,0,0 
  write(21,*) ((fs(i,j,1),i=1,nx,2),j=1,ny,2)
  write(21,*) ((fs(i,j,3),i=1,nx,2),j=1,ny,2)
  write(21,*) ((p(i,j),i=1,nx,2),j=1,ny,2)
  write(21,*) ((T(i,j),i=1,nx,2),j=1,ny,2)
  write(21,*) ((omega(i,j),i=1,nx,2),j=1,ny,2)

  stop 
end program colucci

!*********************************************************************

subroutine rgrid 
  use params
  use coord
  use trnsfrm
  implicit none

  integer :: i, j
  real    :: xmax, ymax, hy, beta, aa, dx, eta(ny)

  dxi=1./(nx-1)
  deta=1./(ny-1)

  !     read externally generated grid 
  !      do 10 j=1,ny 
  !         do 10 i=1,nx        
  !            read(33,*) x(i,j),y(i,j) 
  ! 10   continue

  xmax=120.
  hy=60.  
  ymax=hy/2.
  beta=5. 
  aa=.5/beta*log((1+(exp(beta)-1.)/2.)/(1+(exp(-beta)-1)/2.))
  dx=xmax/(nx-1)
  do j=1,ny 
     do i=1,nx
        eta(j)=(j-1)*deta 
        x(i,j)=(i-1)*dx
        y(i,j)=ymax*sinh(beta*(eta(j)-aa))/sinh(beta*aa)
     end do
  end do


  return
end subroutine rgrid

!*********************************************************************

subroutine metric 
  use params
  use coord
  use trnsfrm
  use delta
  implicit none

  integer :: i, j
  real    :: rjac

  !     inverse metrics 
  !     rmeti(i,j,1)=d(y)/d(eta) 
  !     rmeti(i,j,2)=d(y)/d(xi) 
  !     rmeti(i,j,3)=d(x)/d(eta) 
  !     rmeti(i,j,4)=d(x)/d(xi)          

  !     metrics 
  !     rmet(i,j,1)=d(xi)/d(x) 
  !     rmet(i,j,2)=d(eta)/d(x) 
  !     rmet(i,j,3)=d(xi)/d(y) 
  !     rmet(i,j,4)=d(eta)/d(y) 

  !     compute inverse metrics using 2nd order central difference  
  call ydc2(y(1,1),rmeti(1,1,1)) 
  call xdc2(y(1,1),rmeti(1,1,2)) 
  call ydc2(x(1,1),rmeti(1,1,3)) 
  call xdc2(x(1,1),rmeti(1,1,4)) 

  do j=1,ny 
     do i=1,nx 
        rmeti(i,j,1)=rmeti(i,j,1)/deta  
        rmeti(i,j,2)=rmeti(i,j,2)/dxi  
        rmeti(i,j,3)=rmeti(i,j,3)/deta  
        rmeti(i,j,4)=rmeti(i,j,4)/dxi
     end do
  end do

  !     calculate inverse Jacobian matrix and metrics 
  do j=1,ny 
     do i=1,nx 
        rjaci(i,j)=rmeti(i,j,4)*rmeti(i,j,1) &
             -rmeti(i,j,3)*rmeti(i,j,2)    
        rjac=1./rjaci(i,j) 
        rmet(i,j,1)=rmeti(i,j,1)*rjac  
        rmet(i,j,2)=-rmeti(i,j,2)*rjac  
        rmet(i,j,3)=-rmeti(i,j,3)*rjac  
        rmet(i,j,4)=rmeti(i,j,4)*rjac  
     end do
  end do

  !     compute dx and dy (grid increments) 
  do j=1,ny 
     do i=1,nx-1 
        dx(i,j)=x(i+1,j)-x(i,j)
     end do
  end do
  do j=1,ny 
     dx(nx,j)=dx(nx-1,j)           
  end do

  do j=1,ny-1 
     do i=1,nx 
        dy(i,j)=y(i,j+1)-y(i,j) 
     end do
  end do
  do i=1,nx 
     dy(i,ny)=dy(i,ny-1) 
  end do

  return 
end subroutine metric

!*********************************************************************

subroutine init
  use params
  use coord
  use trnsfrm
  use var
  use strs
  use species
  use vec
  use steps
  use props
  use react
  use same
  implicit none

  integer :: i, j
  real    :: dd, uinfbot, uinftop

  uinftop=1. 
  uinfbot=.25 
  dd=1. 
  do j=1,ny
     do i=1,nx
        u(i,j)=(uinftop+uinfbot)/2. &
             + (uinftop-uinfbot)/2.*tanh(dd*y(i,j))
        v(i,j)=0.0 
        p(i,j)=1./(gamma*M1**2)  
        fs(i,j,1)=.5*(1.+tanh(dd*y(i,j)))
        fs(i,j,2)=1.-fs(i,j,1) 
        fs(i,j,3)=2.*amin1(fs(i,j,1),fs(i,j,2)) 
        if (fs(i,j,1).ge.0.5) then 
           fs(i,j,1)=1.-fs(i,j,3) 
           fs(i,j,2)=0. 
        else 
           fs(i,j,1)=0. 
           fs(i,j,2)=1.-fs(i,j,3) 
        endif
        T(i,j)=1.+Ce*fs(i,j,3) 
        E(i,j)=(T(i,j)-gamma*Ce*fs(i,j,3))/ &
             (gamma*(gamma-1.)*M1**2) + (u(i,j)**2 + v(i,j)**2)/2.
        rho(i,j)=1./T(i,j)  
     end do
  end do

  return
end subroutine init

!*********************************************************************

subroutine advance 
  use params
  use coord
  use trnsfrm
  use var
  use strs
  use species
  use vec
  use steps
  use props
  use react
  use same
  implicit none

  integer :: i, j, k, nstep

  call stress 

  !     predictor step 
  nstep=1

  !     continuity
  do j=1,ny
     do i=1,nx
        A(i,j,1)=rho(i,j)*rjaci(i,j)
        F(i,j)=rho(i,j)*u(i,j)
        G(i,j)=rho(i,j)*v(i,j)
        H(i,j)=0. 
     end do
  end do
  call integ(1,nstep) 

  !     x-momentum
  do j=1,ny
     do i=1,nx
        A(i,j,2)=rho(i,j)*u(i,j)*rjaci(i,j)
        F(i,j)=rho(i,j)*u(i,j)**2 - Txx(i,j)
        G(i,j)=rho(i,j)*u(i,j)*v(i,j) - Txy(i,j)
        H(i,j)=0. 
     end do
  end do
  call integ(2,nstep)

  !     y-momentum
  do j=1,ny
     do i=1,nx
        A(i,j,3)=rho(i,j)*v(i,j)*rjaci(i,j)
        F(i,j)=rho(i,j)*u(i,j)*v(i,j) - Txy(i,j)
        G(i,j)=rho(i,j)*v(i,j)**2 - Tyy(i,j)
        H(i,j)=0. 
     end do
  end do
  call integ(3,nstep)

  !     Energy (total/unit mass)
  do j=1,ny
     do i=1,nx
        A(i,j,4)=rho(i,j)*E(i,j)*rjaci(i,j)
        F(i,j)=(rho(i,j)*E(i,j)-Txx(i,j))*u(i,j) &
             -Txy(i,j)*v(i,j) + qx(i,j)
        G(i,j)=(rho(i,j)*E(i,j)-Tyy(i,j))*v(i,j) &
             -Txy(i,j)*u(i,j) + qy(i,j) 
        H(i,j)=0.  
     end do
  end do
  call integ(4,nstep)

  !     Species equations
  do k=1,nse   
     do j=1,ny
        do i=1,nx
           A(i,j,4+k)=rho(i,j)*fs(i,j,k)*rjaci(i,j)
           F(i,j)=rho(i,j)*u(i,j)*fs(i,j,k) + mfluxx(i,j,k)
           G(i,j)=rho(i,j)*v(i,j)*fs(i,j,k) + mfluxy(i,j,k)
           omega(i,j)=Da1*rho(i,j)**2*fs(i,j,1)*fs(i,j,2)* &
                exp(-Ze*(1./T(i,j) - 1.)) 
           H(i,j)=-omega(i,j)*rjaci(i,j)
        end do
     end do
     call integ(4+k,nstep)
  end do

  !     intermediate (star) values at interior points 
  do j=2,ny-1 
     do i=2,nx-1 
        rho(i,j)=As(i,j,1)/rjaci(i,j)
        u(i,j)=As(i,j,2)/As(i,j,1)
        v(i,j)=As(i,j,3)/As(i,j,1)
        E(i,j)=As(i,j,4)/As(i,j,1)
        do k=1,nse     
           fs(i,j,k)=As(i,j,4+k)/As(i,j,1)
           if (fs(i,j,k).lt.0.) fs(i,j,k)=0.
           if (fs(i,j,k).gt.1.) fs(i,j,k)=1.
        end do
        fs(i,j,3)=1.-(fs(i,j,1)+fs(i,j,2)) 
        T(i,j)=gamma*(gamma-1.)*M1**2* &
             (E(i,j)-(u(i,j)**2+v(i,j)**2)/2.) + Ce*gamma*fs(i,j,3)
        p(i,j)=rho(i,j)*T(i,j)/(gamma*M1**2) 
     end do
  end do

  call bc      

  call stress 

  !     corrector step
  nstep=2

  !     continuity
  do j=1,ny
     do i=1,nx
        A(i,j,1)=rho(i,j)*rjaci(i,j)
        F(i,j)=rho(i,j)*u(i,j)
        G(i,j)=rho(i,j)*v(i,j)
        H(i,j)=0
     end do
  end do
  call integ(1,nstep)

  !     x-momentum
  do j=1,ny
     do i=1,nx
        A(i,j,2)=rho(i,j)*u(i,j)*rjaci(i,j)
        F(i,j)=rho(i,j)*u(i,j)**2 - Txx(i,j)
        G(i,j)=rho(i,j)*u(i,j)*v(i,j) - Txy(i,j)
        H(i,j)=0
     end do
  end do
  call integ(2,nstep)

  !     y-momentum
  do j=1,ny
     do i=1,nx
        A(i,j,3)=rho(i,j)*v(i,j)*rjaci(i,j)
        F(i,j)=rho(i,j)*u(i,j)*v(i,j) - Txy(i,j)
        G(i,j)=rho(i,j)*v(i,j)**2 - Tyy(i,j)
        H(i,j)=0
     end do
  end do
  call integ(3,nstep)

  !     Energy (total/unit mass)
  do j=1,ny
     do i=1,nx
        A(i,j,4)=rho(i,j)*E(i,j)*rjaci(i,j)
        F(i,j)=(rho(i,j)*E(i,j)-Txx(i,j))*u(i,j) &
             -Txy(i,j)*v(i,j) + qx(i,j)
        G(i,j)=(rho(i,j)*E(i,j)-Tyy(i,j))*v(i,j) &
             -Txy(i,j)*u(i,j) + qy(i,j)
        H(i,j)=0. 
     end do
  end do
  call integ(4,nstep)

  !    Species
  do k=1,nse  
     do j=1,ny
        do i=1,nx
           A(i,j,4+k)=rho(i,j)*fs(i,j,k)*rjaci(i,j)
           F(i,j)=rho(i,j)*u(i,j)*fs(i,j,k) + mfluxx(i,j,k)
           G(i,j)=rho(i,j)*v(i,j)*fs(i,j,k) + mfluxy(i,j,k)
           omega(i,j)=Da1*rho(i,j)**2*fs(i,j,1)*fs(i,j,2)* &
                exp(-Ze*(1./T(i,j) - 1.))   
           H(i,j)=-omega(i,j)*rjaci(i,j)
        end do
     end do
     call integ(4+k,nstep)
  end do

  !     new values (n+1) at interior points 
  do j=2,ny-1 
     do i=2,nx-1 
        rho(i,j)=A(i,j,1)/rjaci(i,j)
        u(i,j)=A(i,j,2)/A(i,j,1) 
        v(i,j)=A(i,j,3)/A(i,j,1) 
        E(i,j)=A(i,j,4)/A(i,j,1) 
        do k=1,nse 
           fs(i,j,k)=A(i,j,4+k)/A(i,j,1)
           if (fs(i,j,k).lt.0.) fs(i,j,k)=0.
           if (fs(i,j,k).gt.1.) fs(i,j,k)=1.            
        end do
        fs(i,j,3)=1.-(fs(i,j,1)+fs(i,j,2)) 
        T(i,j)=gamma*(gamma-1.)*M1**2* &
             (E(i,j)-(u(i,j)**2+v(i,j)**2)/2.) + Ce*gamma*fs(i,j,3)
        p(i,j)=rho(i,j)*T(i,j)/(gamma*M1**2) 
     end do
  end do

  call bc 

  return 
end subroutine advance

!*********************************************************************

subroutine integ(nvar,nstep)  
  use params
  use coord
  use trnsfrm
  use var
  use strs
  use species
  use vec
  use veccap
  use steps
  use react
  use props
  use same
  implicit none

  integer, intent(in) :: nvar, nstep

  integer :: i, j, num
  real    :: dFc(nx,ny), dGc(nx,ny)

  call fgcap 
  num=mod(iter,4) 

  !     predictor pass    
  if (nstep.eq.1) then  
     if (num.eq.1) then 
        call xdf4(Fc,dFc) 
        call ydf4(Gc,dGc)  
     elseif (num.eq.2) then 
        call xdb4(Fc,dFc) 
        call ydb4(Gc,dGc)  
     elseif (num.eq.3) then 
        call xdf4(Fc,dFc) 
        call ydb4(Gc,dGc) 
     else 
        call xdb4(Fc,dFc) 
        call ydf4(Gc,dGc) 
     endif

     do j=2,ny-1 
        do i=2,nx-1 
           Aold(i,j,nvar)=A(i,j,nvar) 
           As(i,j,nvar)=A(i,j,nvar) &
                - dt*(dFc(i,j)/dxi+dGc(i,j)/deta) +dt*H(i,j)  
        end do
     end do

     !        corrector pass  
  else 
     if (num.eq.1) then
        call xdb4(Fc,dFc)  
        call ydb4(Gc,dGc) 
     elseif (num.eq.2) then
        call xdf4(Fc,dFc)  
        call ydf4(Gc,dGc)  
     elseif (num.eq.3) then
        call xdb4(Fc,dFc) 
        call ydf4(Gc,dGc) 
     else 
        call xdf4(Fc,dFc) 
        call ydb4(Gc,dGc) 
     endif

     do j=2,ny-1 
        do i=2,nx-1  
           A(i,j,nvar)=.5*(Aold(i,j,nvar) + As(i,j,nvar) &
                -dt*(dFc(i,j)/dxi + dGc(i,j)/deta) +dt*H(i,j))   
        end do
     end do
  endif

  return  
end subroutine integ

!*********************************************************************

subroutine fgcap
  use params
  use vec
  use trnsfrm
  use veccap
  implicit none

  integer :: i, j

  do j=1,ny 
     do i=1,nx 
        Fc(i,j)=rmeti(i,j,1)*f(i,j)-rmeti(i,j,3)*g(i,j)
        Gc(i,j)=rmeti(i,j,4)*g(i,j)-rmeti(i,j,2)*f(i,j) 
     end do
  end do
  return 
end subroutine fgcap

!*********************************************************************

subroutine bc 
  use params
  use coord
  use trnsfrm
  use var
  use strs
  use species
  use vec
  use steps
  use props
  use react
  use same
  implicit none

  integer :: i, j, k
  real    :: uinftop, uinfbot, dd, e1, a1, freq, amplit

  !     zero derivative bc at freestreams (2nd order fwd, back diff) 

  do i=2,nx-1 
     rho(i,1)=(4*rho(i,2)-rho(i,3))/3. 
     u(i,1)=(4*u(i,2)-u(i,3))/3. 
     v(i,1)=(4*v(i,2)-v(i,3))/3. 
     p(i,1)=(4*p(i,2)-p(i,3))/3. 
     E(i,1)=(4*E(i,2)-E(i,3))/3. 
     T(i,1)=(4*T(i,2)-T(i,3))/3. 
     do k=1,nse
        fs(i,1,k)=(4*fs(i,2,k)-fs(i,3,k))/3. 
     end do

     rho(i,ny)=(4*rho(i,ny-1)-rho(i,ny-2))/3. 
     u(i,ny)=(4*u(i,ny-1)-u(i,ny-2))/3.   
     v(i,ny)=(4*v(i,ny-1)-v(i,ny-2))/3.
     p(i,ny)=(4*p(i,ny-1)-p(i,ny-2))/3. 
     E(i,ny)=(4*E(i,ny-1)-E(i,ny-2))/3. 
     T(i,ny)=(4*T(i,ny-1)-T(i,ny-2))/3.  
     do k=1,nse
        fs(i,ny,k)=(4*fs(i,ny-1,k)-fs(i,ny-2,k))/3.  
     end do
  end do

  !     inlet, outlet boundary conditions
  uinftop=1. 
  uinfbot=.25 
  dd=1.
  do j=1,ny
     u(1,j)=(uinftop+uinfbot)/2. &
          +   (uinftop-uinfbot)/2.*tanh(dd*y(1,j))
     v(1,j)=0.0 
     p(1,j)=1./(gamma*M1**2) 
     T(1,j)=T(1,j) 
     fs(1,j,1)=fs(1,j,1) 
     fs(1,j,2)=fs(1,j,2) 
     fs(1,j,3)=fs(1,j,3) 
     E(1,j)=(T(1,j)-gamma*Ce*fs(1,j,3))/(gamma*(gamma-1.)*M1**2) &
          + (u(1,j)**2 + v(1,j)**2)/2.
     rho(1,j)=rho(1,j) 

     !        forcing for v
     e1=.025
     a1=.475*dd
     freq=dd*uinfbot
     amplit=1./cosh(a1*y(1,j))
     v(1,j)=e1*(uinftop+uinfbot)/2.*amplit*sin(freq*time)



     !        zero derivative exit condition (1st order)  
     rho(nx,j)=rho(nx-1,j)
     u(nx,j)=u(nx-1,j)
     v(nx,j)=v(nx-1,j)
     p(nx,j)=p(nx-1,j)
     E(nx,j)=E(nx-1,j)
     T(nx,j)=T(nx-1,j)
     do k=1,nse 
        fs(nx,j,k)=fs(nx-1,j,k)
     end do
  end do

  return 
end subroutine bc

!*********************************************************************

subroutine stress 
  use params
  use coord
  use trnsfrm
  use var
  use strs
  use species
  use vec
  use steps
  use props
  use react
  use same
  implicit none

  integer :: i, j, k
  real    :: dfdx, dfdy, dTdx, dTdy, divV
  real    :: dudx, dudy, dvdx, dvdy, dypdx, dypdy
  real, dimension(nx,ny) :: dudxi, dudeta, dvdxi, dvdeta
  real, dimension(nx,ny) :: dTdxi, dTdeta, dfdxi, dfdeta
  real, dimension(nx,ny) :: phi

  call xdc2(u,dudxi)
  call xdc2(v,dvdxi) 
  call ydc2(u,dudeta) 
  call ydc2(v,dvdeta) 

  do j=1,ny 
     do i=1,nx 
        dudxi(i,j)=dudxi(i,j)/dxi  
        dvdxi(i,j)=dvdxi(i,j)/dxi  
        dudeta(i,j)=dudeta(i,j)/deta  
        dvdeta(i,j)=dvdeta(i,j)/deta  
     end do
  end do

  do j=1,ny 
     do i=1,nx 
        dudx=dudxi(i,j)*rmet(i,j,1)+dudeta(i,j)*rmet(i,j,2)
        dvdx=dvdxi(i,j)*rmet(i,j,1)+dvdeta(i,j)*rmet(i,j,2)
        dudy=dudxi(i,j)*rmet(i,j,3)+dudeta(i,j)*rmet(i,j,4)
        dvdy=dvdxi(i,j)*rmet(i,j,3)+dvdeta(i,j)*rmet(i,j,4)
        divV=dudx+dvdy 
        Txx(i,j)=-p(i,j) + 2.*mu/Re1*(dudx-divV/3.) 
        Txy(i,j)=mu/Re1*(dudy+dvdx)
        Tyy(i,j)=-p(i,j) + 2.*mu/Re1*(dvdy-divV/3.)
     end do
  end do

  do k=1,nse 
     do j=1,ny 
        do i=1,nx 
           phi(i,j)=fs(i,j,k) 
        end do
     end do
     call xdc2(phi,dfdxi)
     call ydc2(phi,dfdeta)
     do j=1,ny 
        do i=1,nx  
           dfdxi(i,j)=dfdxi(i,j)/dxi
           dfdeta(i,j)=dfdeta(i,j)/deta 
        end do
     end do
     do j=1,ny 
        do i=1,nx
           dfdx=dfdxi(i,j)*rmet(i,j,1)+dfdeta(i,j)*rmet(i,j,2)   
           dfdy=dfdxi(i,j)*rmet(i,j,3)+dfdeta(i,j)*rmet(i,j,4)   
           mfluxx(i,j,k)=-mu/Sc/Re1*dfdx
           mfluxy(i,j,k)=-mu/Sc/Re1*dfdy
        end do
     end do
  end do

  call xdc2(T,dTdxi) 
  call ydc2(T,dTdeta) 
  do j=1,ny 
     do i=1,nx 
        dTdxi(i,j)=dTdxi(i,j)/dxi 
        dTdeta(i,j)=dTdeta(i,j)/deta 
     end do
  end do
  do j=1,ny 
     do i=1,nx 
        phi(i,j)=fs(i,j,3) 
     end do
  end do
  call xdc2(phi,dfdxi) 
  call ydc2(phi,dfdeta)
  do j=1,ny 
     do i=1,nx 
        dfdxi(i,j)=dfdxi(i,j)/dxi
        dfdeta(i,j)=dfdeta(i,j)/deta
     end do
  end do
  do j=1,ny 
     do i=1,nx 
        dTdx=dTdxi(i,j)*rmet(i,j,1)+dTdeta(i,j)*rmet(i,j,2)   
        dTdy=dTdxi(i,j)*rmet(i,j,3)+dTdeta(i,j)*rmet(i,j,4)   
        dypdx=dfdxi(i,j)*rmet(i,j,1)+dfdeta(i,j)*rmet(i,j,2)   
        dypdy=dfdxi(i,j)*rmet(i,j,3)+dfdeta(i,j)*rmet(i,j,4)   
        qx(i,j)=-mu*(dTdx-Ce/Le*dypdx)/ &
             ((gamma-1)*M1**2*Pr*Re1)   
        qy(i,j)=-mu*(dTdy-Ce/Le*dypdy)/ &
             ((gamma-1)*M1**2*Pr*Re1)   
     end do
  end do
  return 
end subroutine stress

!*********************************************************************

subroutine tstep 
  use params
  use coord
  use trnsfrm
  use var
  use strs
  use species
  use vec
  use steps
  use props
  use react
  use same
  use delta
  implicit none

  integer :: i, j
  real :: dtn, term1, term2, term3

  dt=1.e20  
  do j=1,ny
     do i=1,nx
        term1=abs(u(i,j))/dx(i,j)   
        term2=abs(v(i,j))/dy(i,j) 
        term3=T(i,j)/M1**2*(1./dx(i,j)**2+1./dy(i,j)**2)
        term3=sqrt(term3)
        dtn=stable/(term1+term2+term3)
        if (dtn.lt.dt) dt=dtn 
     end do
  end do
  return
end subroutine tstep

!*********************************************************************

subroutine xdf4(fcn,dfcn) 
  use params
  implicit none

  integer :: i, j
  real    :: fcn(nx,ny), dfcn(nx,ny) 

  !     4th order forward differencing

  do i=1,nx-2
     do j=1,ny   
        dfcn(i,j)=(-7.0*fcn(i,j)+8.0*fcn(i+1,j)-1.0*fcn(i+2,j))/6. 
     end do
  end do
  do j=1,ny
     dfcn(nx-1,j)=fcn(nx,j)-fcn(nx-1,j)
     dfcn(nx,j)=fcn(nx,j)-fcn(nx-1,j)
  end do

  return
end subroutine xdf4

!*********************************************************************

subroutine xdb4(fcn,dfcn) 
  use params
  implicit none

  integer :: i, j
  real    :: fcn(nx,ny), dfcn(nx,ny) 

  !     4th order backward differencing

  do i=3,nx  
     do j=1,ny 
        dfcn(i,j)=(7.0*fcn(i,j)-8.0*fcn(i-1,j)+1.0*fcn(i-2,j))/6. 
     end do
  end do

  do j=1,ny
     dfcn(1,j)=fcn(2,j)-fcn(1,j)
     dfcn(2,j)=fcn(2,j)-fcn(1,j)  
  end do

  return
end subroutine xdb4

!********************************************************************* 

subroutine ydf4(fcn,dfcn) 
  use params
  implicit none

  integer :: i, j
  real    :: fcn(nx,ny), dfcn(nx,ny) 

  !     4th order forward differencing

  do i=1,nx 
     do j=1,ny-2 
        dfcn(i,j)=(-7.0*fcn(i,j)+8.0*fcn(i,j+1)-1.0*fcn(i,j+2))/6. 
     end do
  end do
  do i=1,nx 
     dfcn(i,ny-1)=fcn(i,ny)-fcn(i,ny-1)  
     dfcn(i,ny)=fcn(i,ny)-fcn(i,ny-1) 
  end do

  return 
end subroutine ydf4

!********************************************************************* 

subroutine ydb4(fcn,dfcn) 
  use params
  implicit none

  integer :: i, j
  real    :: fcn(nx,ny), dfcn(nx,ny) 

  !     4th order backward differencing 

  do i=1,nx 
     do j=3,ny 
        dfcn(i,j)=(7.0*fcn(i,j)-8.0*fcn(i,j-1)+1.0*fcn(i,j-2))/6. 
     end do
  end do
  do i=1,nx 
     dfcn(i,1)=fcn(i,2)-fcn(i,1) 
     dfcn(i,2)=fcn(i,2)-fcn(i,1)
  end do

  return 
end subroutine ydb4

!*********************************************************************

subroutine xdf2(fcn,dfcn)
  use params
  implicit none

  integer :: i, j
  real    :: fcn(nx,ny), dfcn(nx,ny)

  !     2nd order forward differencing

  do i=1,nx-1 
     do j=1,ny
        dfcn(i,j)=fcn(i+1,j)-fcn(i,j)
     end do
  end do
  do j=1,ny
     dfcn(nx,j)=fcn(nx,j)-fcn(nx-1,j)
  end do

  return
end subroutine xdf2

!*********************************************************************

subroutine xdb2(fcn,dfcn)   
  use params
  implicit none

  integer :: i, j
  real    :: fcn(nx,ny), dfcn(nx,ny) 

  !     2nd order backward differencing 

  do i=2,nx 
     do j=1,ny 
        dfcn(i,j)=fcn(i,j)-fcn(i-1,j)
     end do
  end do
  do j=1,ny 
     dfcn(1,j)=fcn(2,j) -fcn(1,j)
  end do

  return 
end subroutine xdb2

!*********************************************************************

subroutine ydf2(fcn,dfcn)   
  use params
  implicit none

  integer :: i, j
  real    :: fcn(nx,ny), dfcn(nx,ny) 

  !     2 nd order forward differencing 

  do i=1,nx 
     do j=1,ny-1 
        dfcn(i,j)=fcn(i,j+1)-fcn(i,j)
     end do
  end do

  do i=1,nx 
     dfcn(i,ny)=fcn(i,ny) -fcn(i,ny-1)
  end do

  return 
end subroutine ydf2

!*********************************************************************

subroutine ydb2(fcn,dfcn)   
  use params
  implicit none

  integer :: i, j
  real    :: fcn(nx,ny), dfcn(nx,ny)

  !     2nd order backward differencing

  do i=1,nx
     do j=2,ny
        dfcn(i,j)=fcn(i,j)-fcn(i,j-1)
     end do
  end do
  do i=1,nx
     dfcn(i,1)=fcn(i,2)-fcn(i,1)
  end do

  return
end subroutine ydb2

!*********************************************************************

subroutine xdc2(fcn,dfcn)
  use params
  implicit none

  integer :: i, j
  real    :: fcn(nx,ny), dfcn(nx,ny)

  !     2nd order central differencing

  do i=2,nx-1
     do j=1,ny
        dfcn(i,j)=.5*(fcn(i+1,j)-fcn(i-1,j))
     end do
  end do
  do j=1,ny 
     dfcn(1,j)=fcn(2,j) -fcn(1,j)
     dfcn(nx,j)=fcn(nx,j) -fcn(nx-1,j)
  end do

  return
end subroutine xdc2

!*********************************************************************

subroutine ydc2(fcn,dfcn)  
  use params
  implicit none

  integer :: i, j
  real    :: fcn(nx,ny), dfcn(nx,ny)

  !     2nd order central differencing

  do i=1,nx
     do j=2,ny-1
        dfcn(i,j)=.5*(fcn(i,j+1)-fcn(i,j-1))
     end do
  end do
  do i=1,nx
     dfcn(i,1)=fcn(i,2) -fcn(i,1)
     dfcn(i,ny)=fcn(i,ny) -fcn(i,ny-1)
  end do

  return
end subroutine ydc2